Hierarchical Liouville-space approach to nonequilibrium dynamical properties of 

quantum impurity systems 
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We propose a hierarchical dynamics approach for evaluation of nonequilibrium dynamic response 
properties of quantum impurity systems. It is based on a hierarchical equations of motion formalism, 
in conjunction with a linear response theory established in the hierarchical Liouville space. This 
provides an accurate and universal tool for characterization of arbitrary response and correlation 
functions of local impurities, as well as transport related response properties. The practicality of our 
proposed approach is demonstrated via the evaluation of various dynamical properties of a single- 
impurity Anderson model. These include the impurity spectral density, local charge fluctuation, local 
magnetic susceptibility, and current-voltage admittance, in both equilibrium and nonequilibrium 
situations. The numerical results are considered to achieve the quantitative accuracy, as long as 
they converge with respect to the truncation of the hierarchy. 

PACS numbers: 71.27.+a, 72.15.Qm, 73.63.Kv 



I. INTRODUCTION 

Recent advances in fabrication, manipulation and mea- 
surement of artificial quantum impurity systems such 
as quantum dots have led to a resurgence of interest 
of nanostructures in both experiment and theory. A 
favorable feature of these nanostructures is the out- 
standing tunability of device parameters. Understand- 
ing the dynamical properties of quantum impurity sys- 
tems is of fundamental importance for the development 
of solid-state quantum information processing 1 - - — and 
single-electron devices^ 5 - Moveover, quantum impurity 
models serve as essential theoretical tools, covering a 
broad range of important physical systems. For instance, 
the Hubbard lattice model can be mapped onto the 
Anderson impurity model via a self-consistent dynam- 
ical mean field theory^— Besides the strong electron- 
electron (e-e) interactions, local impurities are also 
subject to interactions with itinerant electrons in sur- 
rounding bulk materials, which serve as the electron 
reservoir as well as thermal bath. The interplay be- 
tween the local e-e interactions and nonlocal transfer 
coupling gives rise to a variety of intriguing phenom- 
ena of prominent many-body nature, such as Kondo 
effect f 9 - - — Mott metal-insulator transition,— - — and high- 
temperature superconductivity J£r— 

Characterizing the system responses to external per- 
turbation of experimental relevance is of fundamental 
significance in understanding the intrinsic properties of 
quantum impurity systems and their potential applica- 
tions. For instance, the magnetic susceptibility of an 
impurity system reflects the redistribution of electron 
spin under an applied magnetic field, and its investiga- 
tion may have important implications for fields such as 
spintronics. 

For the accurate characterization of dynamical prop- 
erties of the impurity such as the impurity spectral 
function and dynamical charge/magnetic susceptibility, 
a variety of nonperturbative numerical approaches have 



been developed, such as numerical renormalization group 
method ) 11 ' 18 ' 19 density matrix renormalization group 
approach^ - — and quantum Monte Carlo method^ - — 
While most of work has focused on equilibrium proper- 
ties, the accurate characterization of nonequilibrium dy- 
namical properties has remained very challenging. 

In many experimental setups) 9 -^ artificial quantum 
impurity systems attached to electron reservoirs are sub- 
ject to applied bias voltages. This stimulates the ex- 
perimental and theoretical exploration of nonequilibrium 
processes in quantum impurity systems. A variety of in- 
teresting physical phenomena have been observed, which 
originate from the interplay between strong electron cor- 
relation and nonequilibrium dissipation^ - — 

In the past few years, a number of nonper- 
turbative theoretical approaches have been devised 
to treat systems away from equilibrium. These 
include the time-dependent numerical renormaliza- 
tion group methodj^I - — time-dependent density ma- 
trix renormalization group method^ nonequilibrium 
functional renormalization group , 35 ' 36 quantum Monte 
Carlo method*^ - — iterative real-time path integral 
approach ) 40 ' 41 and nonequilibrium Bethe ansatz42r— 
Despite the progress made, quantitative accuracy is not 
guaranteed for the resulted nonequilibrium properties, 
because of the various simplifications and approximations 
involved in these approaches. Therefore, an accurate 
and universal approach which is capable of addressing 
nonequilibrium situations is highly desirable. 

In this work we propose a hierarchical dynamics ap- 
proach for the characterization of nonequilibrium re- 
sponse of local impurities to external fields. A gen- 
eral hierarchical equations of motion (HEOM) approach 
has been developed^ 5 - - — which describes the reduced 
dynamics of open dissipativc systems under arbitrary 
time-dependent external fields. The HEOM theory re- 
solves the combined effects of e-e interactions, impurity- 
reservoir dissipation, and non-Markovian memory in a 
nonperturbative manner. In the framework of HEOM, 
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the nonequilibrium dynamics are treated by following the 
same numerical procedures as in equilibrium situations. 
The HEOM theory is in principle exact for an arbitrary 
equilibrium or nonequilibrium system, provided that the 
full hierarchy inclusive of infinite levels are taken into 
account.— In practice, the hierarchy needs to be trun- 
cated at a finite level for numerical tractability. The 
convergence of calculation results with respect to the 
truncation level should be carefully examined. Once the 
convergence is achieved, the numerical outcome is consid- 
ered to reach quantitative accuracy for systems in both 
equilibrium and nonequilibrium situations. 

It has been demonstrated that the HEOM approach 
leads to an accurate and universal characterization of 
strong electron correlation effects in quantum impurity 
systems, and treats the equilibrium and nonequilibrium 
scenarios in a unified manner. For the equilibrium prop- 
erties of Anderson model systems, the HEOM approach 
achieves the same level of accuracy as the latest state- 
of-the-art NRG method^ In particular, the universal 
Kondo scaling of zero-bias conductance and the logarith- 
mic Kondo spectral tail have been reproduced quantita- 
tively. For systems out of equilibrium, numerical calcula- 
tions achieving quantitative accuracy remain very scarce. 
One of the rare cases where numerically exact solution 
is available is the dynamic current response of a nonin- 
teracting quantum dot to a step-pulse voltagei 49 ' 50 This 
has been precisely reproduced by the HEOM approach^ 
However, there are very few calculations at the level of 
quantitative accuracy for systems involving strong e-e in- 
teractions, since most of the existing methods involve in- 
trinsic approximations. Based on the HEOM formalism, 
quantitative accuracy should be achieved once the nu- 
merical convergence with respect to the truncation level 
of hierarchy is reached. 

There are two schemes to evaluate the response 
properties of quantum impurity systems in the frame- 
work of HEOM: (i) Calculate relevant system correla- 
tion/response functions based on a linear response the- 
ory constructed in the HEOM Liouvillc spaced and (ii) 
solve the EOM for a hierarchical set of density opera- 
tors to obtain the transient reduced dynamics of system 
in response to time-dependent external perturbation, fol- 
lowed by a finite difference analysis. These two schemes 
are completely equivalent in the linear response regime, 
as have been verified numerically. In previous studies, 
we had employed the above second scheme to evaluate 
the dynamic admittance (frequency-dependent electric 
current in response to external voltage applied to cou- 
pling electron reservoirs) of quantum dot systems, which 
had led to the identification of several interesting phe- 
nomena, including dynamic Coulomb blockado^ and dy- 
namic Kondo transition^ and photon-phonon-assisted 
transport^ 

In this work, we will elaborate the above first scheme 
of HEOM approach. The external perturbation may as- 
sociate with an arbitrary operator in the impurities sub- 
space, or originates from a homogeneous shift of elec- 
trostatic potential (and hence the chemical potential) of 
electron reservoir. The detailed numerical procedures 
will be exemplified through the evaluation of a variety of 
response properties of a single-impurity Anderson model, 



including the impurity spectral density function, local 
charge fluctuation spectrum, local magnetic susceptibil- 
ity, and dynamic admittance. 

The remainder of paper is organized as follows. We will 
first give a brief introduction on the HEOM method in 
Sec. mi In Sec. lIIII we will elaborate the establishment of a 
linear response theory in the HEOM Liouvillc space. Cal- 
culation on system correlation/response functions which 
are directly relevant to the response properties of primary 
interest will be discussed in detail. We will then provide 
numerical demonstrations for the evaluation of various 
dynamical properties in Sec. lIVI Finally, the concluding 
remarks will be given in Sec.fVl 

II. A REAL-TIME DYNAMICS THEORY FOR 
NONEQUILIBRIUM IMPURITY SYSTEMS 

A. Prelude 

Consider a quantum impurity system in contact with 
two electron reservoirs, denoted as the a = L and R 
reservoirs, under the bias voltage V = /iL — /Mr- The 
total Hamiltonian of the composite system assumes the 
form of 

-fftotal = H sys + ^(e Q fe + fl a ) d^dak 
a k 

+ ]T (tc^dl^ + h.c.) . (i) 

The impurity system Hamiltonian H sys is rather gen- 
eral, including many-particle interactions and external 
field coupling. Its second quantization form is given in 
terms of electron creation and annihilation operators, 
aj, = at and a„ = a^, which are associated with the 
system spin-state [i. The reservoirs are modeled by a 
noninteracting Hamiltonian; see the second term on the 
right-hand side (rhs) of Eq. (fT|), where d} ak (d a k) and e a k 
are the creation (annihilation) operator and energy of 
single-electron state |fc) electron of a- reservoir, respec- 
tively. While the equilibrium chemical potential of to- 
tal system is set to be = 0, the reservoir states 
are subject to a homogeneous shift, under applied 
voltages. The last term on the rhs of Eq. ([1]) is in 
a standard transfer coupling form, which is responsi- 
ble for the dissipative interactions between the system 
and itinerary electrons of reservoirs. It can be recast as 

h' = E QM (/i«M + & i whcrc /i = = 

(/~ M ) f . Throughout this paper we adopt the atomic unit 
e = h = 1 and denote /3 = 1/(/cbT), with ks being the 
Boltzmann constant and T the temperature of electron 
reservoirs. Introduce also the sign variables, a = +/ — 
and a = —a the opposite sign of a. 

The a-reservoir is characterized by the spectral den- 
sity J aflv (u) = nJ2k t ^kf 1 , t <xki>S(u - e ak ). It influences 
the dynamics of reduced system through the reservoir 
correlation functions C^f v {t - r) = {/^(i)/^(r))«, 
Here, ((-)) a ee tr„ [(•) e -^-]/tr a ( e -^-) and f^(t) ^ 
e tHat f^ fl e~ at , with H a being the Hamiltonian of a- 
reservoir. The superscript "st" highlights the stationary 
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feature of the nonequilibrium correlation function, under 
a constant p a . It is related to the reservoir spectral den- 
sity, J afJ ,„(uj) = J~^(uj) = J£ vll {u), via the fluctuation- 
dissipation theorem^ 5 - 



duj- 



1 



(2) 



Physically, C£jf*(t), with a = + or — , describes the pro- 
cesses of electron tunneling from the a-reservoir into the 
specified system coherent state or the reverse events, re- 
spectively. 

We will be interested in nonequilibrium dynamic re- 
sponses to a time-dependent external field acting on ei- 
ther the local system or the reservoirs. For the latter case, 
we include a time-dependent shift in chemical potential 
SA a (t), on top of the constant p a , to the a-reservoir. Its 
effect can be described by rigid homogeneous shifts for 
the reservoir conduction bands, resulting in the nonsta- 
tionary reservoir correlation functions of 



dt' 5A a (t') 



(3) 



This is the generalization of C%*(t) = e <TV °'(7^(i), 
as inferred from Eq. ([2]), with the equilibrium counter- 
part being of /i° q = 0. In the following, we focus 
on the situation of diagonal reservoir correlation, i.e., 
= ^»<W> and C£#(t) = C^(t)S^. In 
constructing closed HEOM^ we should expand C^At) 
in a finite exponential series, 



(t) 



M 
rn — 1 



(4) 



Involved are a total number of M = N' + N poles from 
the reservoir spectral density and the Fermi function 
in the contour integration evaluation of Eq. Vari- 
ous sum-over-poles schemes have been developed, includ- 
ing the Matsubara spectrum decomposition scheme*^ 5 - 
a hybrid spectrum decomposition and frequency dis- 
persion scheme^ the partial fractional decomposition 
scheme^ and the Pade spectrum decomposition (PSD) 
scheme ) 55 ' 56 with the primary focus on the Fermi func- 
tion. To our knowledge, the PSD scheme has the best 
performance until now. We will come back to this is- 
sue later; see the remark- (6) in Sec. Ill Bl In the present 
work we use the [N—l/N] PSD scheme.—^ 6 - It leads to a 
minimum M = N' + N in the exponential expansion of 
Eq. (fj| and thus an optimal HEOM construction^ 5 - 

The exponential expansion form of the reservoir corre- 
lation function in Eq. ((4]) dictates the explicit expressions 
for the HEOM formalism^ 5 - For bookkeeping we intro- 
duce the abbreviated index j = {aapm} for jj = Ja^m 
and so on, or j = {up} for aj = a° Denote also 
j = {ffa/im} or {ap} whenever appropriate, with a be- 
ing the opposite sign of a = + or — . The dynamical 
variables in HEOM are a set of auxiliary density opera- 
tors (ADOs), {Pj?.. jn (t);n = 0,1,- •■ ,L}, with L being 
the terminal or truncated tier of hierarchy. The zeroth- 
tier ADO is set to be the reduced system density ma- 
trix, p (0) (i) = p(t) = tr ba th [ptotal(*)], i-e-, the trace of 



the total system and bath composite density matrix over 
reservoir bath degrees of freedom. 



B. Hierarchical equations of motion formalism 

The HEOM formalism has been constructed from the 
Feynman- Vernon influence functional path integral the- 
ory, together with the Grassmann algebra^ 5 - The initial 
system-bath decoupling used for expressing explicitly the 
influence functional is set at the infinite past. It does not 
introduce any approximation for the characterization of 
any realistic physical process starting from a stationary 
state, which is defined via the HEOM that includes the 
coherence between the system and bath. The detailed 
construction of HEOM is referred to Ref. |45|. Here we 
just briefly introduce the HEOM formalism and discuss 
some of its key features. 

The final HEOM formalism readsi 5 - 



.(«) 

Ph'-jn 



(5) 



The boundary conditions are = p^~^ = 0, together 
with a truncation by setting all p( Tl>L ' = 0. The ini- 
tial conditions to Eq. ([S]) will be specified in conjunction 
with the evaluation of various response and correlation 
functions in Sec. Mil 

The time-dependent damping parameter 7j"...j n (i) hi 
Eq. (|5|) collects the exponents of nonstationary reservoir 
correlation function [cf. Eqs. ^ and @]: 



Eb,.-^A„(t)] a 



(6) 



This expression has been used directly in the HEOM 
evaluation of transient current dynamical properties un- 
der the influence of arbitrary time-dependent chemical 
potentials applied to electrode leads i 46 i 51 ' 52 i 57 Note that 
7j = 7a M m = 1%$m- In Sec.|HICJ we will treat 

5A a (t) as perturbation and derive the corresponding lin- 
ear response theory formulations for various transport 
current related properties under nonequilibrium (p a 0) 
conditions. 

To evaluate nonequilibrium correlation functions of lo- 
cal system via linear response theory (cf. Sec. lIII Bl . the 
time-dependent reduced system Liouvillian in Eq. ([5J is 
assumed formally the form of 



C(t) = £ s + 6C(t). 



(7) 



Here, C s ■ = [H sys , •] remains the commutator form in- 
volving two -f/gys-actions onto the bra and ket sides in- 
dividually. However, the time-dependent perturbation 
5C(t) may act only on one side, in line with the HEOM 
expressions for local system correlation functions 1 48 i 58 ~ ^l 
Other features of HEOM and remarks, covering both 
the theoretical formulation and numerical implementa- 
tion aspects, are summarized as follows. 
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(1) The Fermi- Grassmannian properties: (i) All j- 

indexes in a nonzero n* h -tier ADO, pj"l.j , must be 
distinct. Swap in any two of them leads to a minus 

(2) (2) 

sign, such as p) = —p) . In line with this property, 
the sum of the tier-up dependence in Eq. ([5]) runs only 
over j ^ {ji,'" ,jn\'i (H) Involved in Eq. J5]) are also 
Aj = AZ and Cj = C* m - They are Grassmann superop- 
erators, defined via their actions on an arbitrary operator 
of fermionic or bosonic (bi-fcrmion) nature, O f or O b , by 

^,O f/b = a ? 6 F/B T O p/B a ? , 

J . J „ / (8) 

CjO F/B = %%O p/B ± VjO F/B a 3 ■ 

In particular, even-tier ADOs are bosonic, while odd-tier 
ones are fermionic. The case of opposite parity would also 
appear in conjunction with applications; see comments 
following Eq. ([15]). 

(2) Physical contents of ADOs: While the zero-tier 
ADO is the reduced density matrix, i.e., p'°'(t) = p(t), 

the first-tier ADOs, p^ = p a alim , are related to the elec- 
tric current through the interface between the system and 
a-reservoir, I a (t), as follows, 

I a (t) = -2 Im Y, Tr Hp^ m (*)] • (9) 

fxm 

Moreover, we have £ m p£ Mm (i) = tr ba th[/S M (>)ptotal(i)], 
and can further relate trbath[/a^(*)/av(*)Ptotai(t)] to 
the second-tier ADOs, and so on. Note that = 

^'(EiWi)^"' = [f*u(t)] 1 are defined in the 
bath-space only. Apparently, the Fermi- Grassmannian 
properties in remark-(l) above are rooted at the 
fermionic nature of individual {fa^}- 

(3) Hcrmitian property: The ADOs satisfy the Hcrmi- 

tian relation of [pj"... 3 - n (*)] = Pj (*)) whenever the 
perturbed idC(t) action assumes Hcrmitian; see the com- 
ments following Eq. (JT)). 

(4) Nonperturbative nature: The HEOM construction 
treats properly the combined effects of system-bath cou- 
pling strength, Coulomb interaction, and bath memory 
time scales, as inferred from the following observations. 
(i) For noninteracting electronic systems, the coupling 
hierarchy stops at second tier level (L = 2) without 
approximation^ 5 - (ii) HEOM is of finite support, con- 
taining in general only a finite number of ADOs. Let 
K be the number of all distinct j-indexes. Such a num- 
ber draws the maximum tier level L max = K, at which 
the HEOM formalism ultimately terminates. The to- 
tal number of ADOs, up to the truncated tier level L, 
is J2t=o n \(K-n)\ ^ as L < L max = K; (iii) The 
hierarchy resolves collectively the memory contents, as 
decomposed in the exponential expansion of bath corre- 
lation functions of Eq. Q . It goes with the observation 

that an individual ADO, Pj\ ..j n , is associated with the 
collective damping constant Re 7^™. • in Eq. ([5]). Mean- 
while Pj™..j has the leading (2n) th -order in the overall 
system-bath coupling strength. One may define proper 
non-Marvokianicity parameters to determine in advance 



the numerical importance of individual ADOsj£i~— (iv) 
Convcrgcncy tests by far - For quantum impurity sys- 
tems with nonzero e-e interactions, calculations often 
converge rapidly and uniformly with the increasing trun- 
cation level L. Quantitatively accurate results are usually 
achieved at a relatively low value of L. 

(5) Nonequilibrium versus equilibrium: The HEOM 
formalism presented earlier provides a unified approach 
to equilibrium, nonequilibrium, time-dependent and 
time-independent situations. In general, the number K 
of distinct ADO indexes amounts to K = 2N a N IJj M , as 
inferred from Eq. (j4|), with being the number spin- 
orbitals of system in direct contact to leads. The factor 
2 accounts for the two choices of the sign variable a, while 
N a = 2 for the distinct a = L and R leads. Interestingly, 
in the equilibrium case, together with the Jl(w) oc Jr(w) 
condition, one can merger all leads into a single lead to 
have the reduced K = 2N fl M. The resulting equilibrium 
HEOM formalism that contains no longer the a-index 
can therefore be evaluated at the considerably reduced 
computational cost. 

(6) Control of accuracy and efficiency: The bath cor- 
relation function in exponential expansion of Eq. (j4]) dic- 
tates the accuracy and efficiency of the HEOM approach, 
(i) The accuracy in the exponential expansion of Eq. (2J) 
is found to be directly transferable to the accuracy of 
HEOM. In other words, HEOM is exact as long as the 
expansion is exact; (ii) The expansion of Eq. ^ is uni- 
formly convergent, and becomes exact when M goes to 
infinity, for any realistic bath spectral density with fi- 
nite bandwidth at finite temperature (T ^ 0); (iii) The 
[N—1 /N] PSD scheme adopted in this work is considered 
to be the best among all possible sum-over-poles expan- 
sion of Fermi function! 55 ' 56 ' 64 In particular it is dramati- 
cally superior over the commonly used Matsubara expan- 
sion expression. The PSD scheme leads to the optimal 
HEOM, with a minimum K-space size, for either equilib- 
rium or nonequilibrium discussed in remark- (5) 
above. 

(7) Computational cost: The CPU time and memory 
space required for HEOM calculations are rather insen- 
sitive to the Coulomb coupling strength and to the equi- 
librium versus nonequilibrium and time-dependent and 
time- independent types of evaluations. However, it grows 
exponentially as the temperature T —y 0, with respect 
to system-bath hybridization strength, due to the signifi- 
cant increase of both the converged if-space and L-space 
sizes. 

To conclude, HEOM is an accurate and versatile tool, 
capable of universal characterizations of real-time dy- 
namics in quantum impurity systems, in both equilib- 
rium and nonequilibrium cases. These remarkable fea- 
tures have been demonstrated recently in several com- 
plex quantum impurity systems^ with the focus mainly 
on equilibrium properties. The HEOM approach is also 
very efficient. Calculations often converge rapidly and 
uniformly with the increasing truncation level L. Quan- 
titatively accurate results are usually achieved at a rela- 
tively low level of truncation^ We will show in Scc. lIVI 
that these features will largely remain in the evaluations 
of nonequilibrium properties. 
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III. NONEQUILIBRIUM RESPONSE THEORY 
A. Linearity of the hierarchical Liouville space 

To highlight the linearity of HEOM, we arrange the 
involving ADOs in a column vector, denoted symbolically 

as 

p(t) = {p(t),pfXt),p?l(t),---}. (10) 

Thus, Eq. ([5]) can be recast in a matrix- vector form (each 
element of the vector in Eq. (fT0| is a matrix) as follows, 

p = -i£(t)p, (11) 

with the time-dependent hierarchical-space Liouvillian, 
as inferred from Eqs. ©-([7]), being of 

C(t) = C s + 6C(t)X + SV(t) . (12) 

It consists not just the time-independent C s part, but 
also two time-dependent parts and each of them will be 
treated as perturbation at the linear response level soon. 
Specifically, 5C{t)X, with X denoting the unit operator in 
the hierarchical Liouville space, is attributed to a time- 
dependent external field acting on the reduced system, 
while SV(t) is diagonal and due to the time-dependent 
potentials SA a (t) applied to electrodes. 

We may denote 5A a (t) = x a SA(t), with < xl = 
1 + x R < 1; thus 8A(t) = 8A L (t) - SA R (t). It specifies 
the additional time-dependent bias voltage, on top of the 
constant V = /L*l — /^r, applied across the two reservoirs. 
As inferred from Eq. ([5]) , we have then 

6V(t) = -S6A(t), (13) 
where S = diag{0, S^)„j n ;n = 1, • • ■ , L}, with 

n 
r=l 

Note that S m = 0. 

The additivity of Eq. (JT^J) and the linearity of HEOM 
lead readily to the interaction picture of the HEOM dy- 
namics in response to the time-dependent external distur- 
bance SC(t) = SC(t)X + 5V(t). The initial unperturbed 
ADOs vector assumes the nonequilibrium steady-state 
form of 

p»\T,V) = {p,pf\pfl,---}, (15) 

under given temperature T and constant bias voltage V. 
It is obtained as the solutions to the linear equations, 
C s p st (T, V) = 0, subject to the normalization condition 
for the reduced density matrix i 45 i 46 i 51 The unperturbed 
HEOM propagator reads Q s {t) = cxp(—iC s t). Based 
on the first-order perturbation theory, Sp(t) = p(t) — 
p st {T 7 V) is then 

8p(t) = -i[ dTg s (t-T)SC(T)p st (T,V). (16) 
Jo 

The response magnitude of a local system observ- 
able A is evaluated by the variation in its expectation 



value, 5A(t) = Ti{A5p(t)}. Apparently, this involves the 
zcroth-ticr ADO 5p{t) in Sp(t) = {8p(t), Sp^}.. jn (t); n = 
1, ■ • • , £}. In contrast, the response current under ap- 
plied voltages cannot be extracted from 6p(t), because 
the current operator is not a local system observable. 
Instead, as inferred from Eq. while the steady-state 
current I a through a-reservoir is related to the steady- 
state first-tier ADOs, pj = Pa^rm the response time- 
dependent current, 5I a (t) = I a (t) — I a , is obtained via 
5pf\t) = Sp^ m (t) : 

The above two situations will be treated respectively, 
by considering 5C(t) = 8C(t)X and 8C(t) = SV(t), in 
the coming two subsections: Scc. lIIIBl treats the local 
system response to a time-dependent external field act- 
ing on the reduced system, while Sec. lIII CI addresses the 
issue of electric current response to external voltage ap- 
plied to reservoirs. 



B. Nonequilibrium correlation and response 
functions of system 

Let A and B be two arbitrary local system operators, 
and consider the correlation functions, Cab{J ~ T ) = 
(A(t)B( T )) st and S AB (t - r ) = ({A(t), B(r)}) st , and re- 
sponse function, XABit — r) = i([A(t), B(r)]) st . It is well 
known that for the equilibrium case they are related to 
each other via the fluctuation-dissipation theorem. The 
nonequilibrium case is rather complicated, and the re- 
lation between nonequilibrium correlation and response 
functions is beyond the scope of the present paper. 

We now focus on the evaluation of local system cor- 
relation/response functions with the HEOM approach. 
This is based on the equivalence between the HEOM- 
space linear response theory of Eq. (fT6"| and that of the 
full systcm-plus-bath composite space. 

We start with the evaluation of nonequilibrium steady- 
state correlation function Cab (t) = (A(t)B(0)) st , as fol- 
lows. By definition, the system correlation function can 
be recast into the form of 

C AB {t) = TT total {Ag to Ut)[Bpf otal (T,V)}} 
= Tr tota i[A5totai(i)] 

= Tr[Ap(t)]. (17) 

The Ptotai^'^O an d £7totai(i) in the nrs t identity arc 
the steady-state density operator and the propagator, 
respectively, in the total system-bath composite space 
under constant bias voltage V. Define in the last two 
identities of Eq. ([T7J are also Ptotai(i) = £total(*)p t otal(0) 
and p{t) = tr b athPtotal(*), with /5 to tai(0) = Bpf otal (T,V). 
Equation (|1T[) can be considered in terms of the linear 
response theory, in which the perturbation Liouvillian 
induced by an external field 8e(t) assumes the form of 
— i8£(t)(-) = B(-)5e(t), followed by the observation on 
the local system dynamical variable A. Both A and B 
can be non-Hcrmitian. Moreover, 5C{t) is treated for- 
mally as a perturbation and can be a one-side action 
rather than having a commutator form. 
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For the evaluation of CUb(£) with the HEOM-space 
dynamics, the corresponding perturbation Liouvillian is 
SC(t) = 5C(t)X, with the above defined 5C(t). It 
leads to -i5C{T)p st (T,V) = Bp st {T,V)6e(T) involved 
in Eq. (|T5)) . The linear response theory that leads to the 
last identity of Eq. (fT7)) is now of the p(t) being just the 
zeroth-tier component of 

p(t) = {p(t), pf\t), pf] 2 (t), ■■■} = a (t)p(O), (18) 

with the initial value of [cf. Eq. t|15[t] 

p(0) = Bp st (T, V) = {B Pl Bpf\ Bpfl ,.-■}. (19) 

The HEOM evaluations of Sab (t) and xab(J) are similar, 
but with the initial ADOs of p(0) = {B, p st (T, V)} and 
i[B, p st (T, V)], respectively. 

Care must be taken when propagating Eq. (fT8|) . for the 
HEOM propagator Q s (t) involving the Grassmann su- 
peroperators Aj and Cj defined in Eq. (|8) . Note that 
the steady-state system density operator p is always of 
the Grassmann-even (or bosonic) parity. Therefore, the 
zeroth-tier ADO pit) in the above cases is of the same 
Grassmann parity as the operator B, while the ADOs at 
the adjacent neighboring tier level are of opposite par- 
ity. The HEOM propagation in Eq. ([18} is then specified 
accordingly. 

It is also worth pointing out that the HEOM evaluation 
of equilibrium correlation and response functions of the 
local system can be simplified when Jl(w) oc Jr(w). In 
this case, two reservoirs can be combined as a whole en- 
tity bath, resulting in a HEOM formalism that depends 
no longer on the reservoir-index a. 

C. Current response to applied bias voltages 

1. Dynamic differential admittance 

Consider first the differential circuit current through 
a two-terminal transport system composed of an quan- 
tum impurity and two leads, 5I(t) = ^[SIh(t) — SIs.(t)], 
in response to a perturbative shift of reservoir chemical 
potential 5A(t). 

We have 

6I a (t) = [ dr G a (t — r) SA(t). (20) 
Jo 

The HEOM-space dynamics results in 

G a (t)=2Re^Tr[a+p- Mm (t)], (21) 

ft in 

with Paum(t) denoting the first-tier ADOs in p{t) 
[Eq. dTSJ)J with the initial value of [cf. Eqs. (p^]) - (fT5|) ] 

p(0) = ~Sp s \T,V) ee -{0, S^pf\sflpfl,- ■ ■}. (22) 
Denote the half-Fourier transform, 

G a {oj) ee / dte lult G a (t). (23) 
Jo 

The admittance is given by G(lo) = ^[G^u) — Gr(w)], 
with its zero- frequency component recovering the steady- 
state differential conductance as dl/dV = G(ui = 0). 



2. Current-number and current- current response functions 

Consider now the differential current SI a (t) in re- 
sponse to an additional time-dependent chemical poten- 
tial SA a '(t) applied on a specified a'-reservoir. Note 
that the Hamiltonian of the total composite system, 
Eq. (Q]), is now subject to a perturbation of 5H to ted(t) = 
N a '5A a >(t), with N a > = J^k^a'k^a'k being electron 
number operator of the a'-rcservoir. Thus, the hierar- 
chical Liouvillc space linear response theory leads to 

SI a (t) = f dr G aa , (t - t) 8A a , (r) , (24) 
Jo 

where the kernel is characterized by the nonequilibrium 
steady-state current-number response function, 

G aal (t-r) = -i{[I a (t),N a ,(T)]) st , (25) 

with ((.)) rt = Tr total [(-)p s 4(T,V)}. Equation 023 can 
be derived by following Eqs. (j9|), (fT6|) and (|24|) . and the 
HEOM evaluation of G aa i (t — r) can be achieved as fol- 
lows. Equation (fT3j) is recast as SC(t) = —S a >SA a i(t), 
where S a > = diag{0, . . - n } , where Sj 1 ...j n is similar to 
S£}.. jn of Eq. (JHJ) but with x a = 5 aa ,. Therefore, 

n 

S h-jn = E( CT ^«')^V • (26) 
r=l 

Its rhs collects the signs (er = +1 or — 1) in the involving 
(j ee {cra/im})-indexes whenever a = a'. The suitable 
initial values for the vector of ADOs are 

MO) = S a ,p st (T, V) = - {0, Sfpf\ S<pfl,- ■ ■} , 

followed by the unperturbed HEOM-space evolution, 

p a >(t) = g s (t)p a ,(0) = {p(t;a'), pf\t;a'),---}. (27) 

The involving first-tier ADOs, p^\t]a') = Pa^m{^ a ')i 
are used to evaluate the current-number response func- 
tion [cf. Eq. pljl]: 

G aa ,(t) = 2Rc^Tr [a+~p^ m (t; a 1 )] . (28) 
fj,m 

Apparently, G a (t) = x^G a i,(t) + x^G a ^i(t), which is just 
the dynamic admittance considered in Sec. lIII C T1 

The nonequilibrium steady-state current-current re- 
sponse function, Xaa' (t) , can be obtained numerically by 
taking the time derivative of G aa '(t), 

Xaa >{t) = i([Ia(t)Ja'(0)]U = G aa ,(t). (29) 

In the hierarchical Liouville space, Xaa' (t) can be explic- 
itly expressed by the zeroth-, first- and second-tier ADOs, 
as inferred from Eq. ([28} and the EOM for Puumfa a ') ■ 
Its Fourier transform, the current-current response spec- 
trum, may carry certain information about the shot noise 
of the impurity system. 

In general, the correlation/response functions between 
an arbitrary local system operator A and the electron 
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number operator N a > of the a'-electrode can be evalu- 
ated via the zeroth-tier ADO p(t;a') of Eg. (|27p . such 
as i{[A(t),N a >(0)]) Bt = -Tr[Ap(t;a% by using the 
HEOM Liouville propagator. Its time derivative gives 

i([A(t)J a ,(o)]) 8t . 



IV. RESULTS AND DISCUSSIONS 

We now demonstrate the numerical performance of 
the HEOM approach on evaluation of nonequilibrium 
response properties of quantum impurity systems. The 
hierarchical Liouville-space linear response theory estab- 
lished in Sec. IIIII is employed to obtain the relevant corre- 
lation/response functions, from which the response prop- 
erties are extracted. 

It is worth emphasizing that the numerical examples 
presented in this section aim at verifying the accuracy 
and universality of the proposed methodology, rather 
than addressing concrete physical problems. To this end, 
the widely studied standard single-impurity Anderson 
model (SIAM) is considered. The Hamiltonian of the 
single-impurity is -ff sys = e-j-n-j- + e^n^ + Un^n^, with 
fin = aj^dfj, being the electron number operator for the 
spin-/! (j" or J,) impurity level. The impurity is coupled 
to two noninteracting electron reservoirs (a = L and R). 
For simplicity, the spectral (or hybridization) function of 
a- reservoir assumes a diagonal and Lorentzian form, i.e., 

J af iv(u) = 6p V / M _^ ff a , with T a and W a being the 
linewidth and bandwidth parameters, respectively. 

Note that the same set of system parameters are 
adopted for all calculations (except for specially spec- 
ified): e = e t = e; = -0.5, U = 1.5, T = 0.02, 
r = T L = T R = 0.1, W h = W K = 2, all in units of 
meV. The nonequilibrium situation concerns a steady 
state defined by a fixed bias voltage applied antisymmet- 
rically to the two reservoirs, i.e., /!l = — /!r = y with 
A = — V = 0.2, and/or 0.7 meV. A recently developed 
[N—l/N] Pade spectrum decomposition schem e 55 ' 56 with 
N = 8 (i.e., M = 9) is used for the efficient construction 
of the hierarchical Liouville propagator associated with 
Eq.©. 

To obtain quantitatively accurate numerical results, we 
increase the truncation level L and the number of expo- 
nential terms M continually until convergence is reached. 
Table U lists the probabilities that the impurity is singly 
occupied by spin-/z electron (P M = (p,\p(T, V)\p) with 
/i =t or |); or doubly occupied (P n = (U\p(T, V)\U))- 
Here, p(T, V) is the nonequilibrium steady-state reduced 
density matrix under temperature T and antisymmet- 
ric applied voltage V. Calculations are done at differ- 
ent truncation level L (up to L = 5) and fixed M = 9. 
Apparently, the HEOM results converge rapidly and uni- 
formly with the increasing L, i.e., with higher-tier ADOs 
included explicitly in Eq. ([5]). In particular, the remain- 
ing relative deviations between the results of L = 4 and 
L = 5 are less than 0.1%, indicating that the L = 4 level 
of truncation is sufficient for the present set of parame- 
ters. It is also affirmed M = 9 is sufficient to yield con- 
vergent results; see Supplemental Material^ 6 - These are 
further affirmed by the calculated steady-state current 



L 




Pn 


I (nA) 


1 


0.500 (0.500) 


0.001 (0.000) 


0.003 


2 


0.441 (0.462) 


0.025 (0.027) 


4.654 


3 


0.439 ( 0.454) 


0.024 (0.025) 


4.920 


4 


0.440 (0.457) 


0.024 (0.024) 


4.799 


5 


0.440 (0.457) 


0.024 (0.024) 


4.799 



TABLE I: Spin-/! single- and double-occupation probabilities 
(P M and -Pfi)i an d steady-state current of an SIAM with two 
electrons reservoirs under an antisymmetrically applied bias 
voltage of A = —V = 0.7 meV. Calculations are done by 
solving the HEOM of Eq. © truncated at different level L. 
The parameters are adopted are (in units of meV): e — — 
£4. = -0.5, U = 1.5, T L = r R = 0.1, W L = Wr = 2, and T = 
0.02. For comparison, the numbers of equilibrium situation 
of V = are shown in the parentheses. 

I(V) across the impurity, which also converges quantita- 
tively with rather minor residual uncertainty at L = 4 
and M = 9. Note also that the truncation at L = 1 level 
results in the sequential current contribution, which is 
negligibly small for the present nonequilibrium system 
setup. The values of / evaluated at different truncation 
levels clearly indicate the current contributions from dif- 
ferent orders of cotunneling processes. 

In the following, we first show the spectral function of 
the SIAM calculated by using the HEOM approach (see 
Fig-UJi an d then present the evaluation of some typi- 
cal response properties in both equilibrium and nonequi- 
librium situations. These will include the local charge 
fluctuation spectrum Sq(lo) (see Fig. [2]), local magnetic 
susceptibility ^mH (see Fig-EJ) , and differential admit- 
tance G(lu) (see Fig.[5|)- All calculations are carried out 
at the truncation level of L = 4 and M = 9. Based on 
the analysis of Table HI the resulting response properties 
are expected to be quantitatively converged with respect 
to L = 4 and M = 9. 

Figure [T] depicts the HEOM calculated spin-/! spectral 
function of the impurity, 

A„M = i Rc [l°°dte iut ({a,(t),al(0)}) st } . (30) 

The effect of bias voltage A on A^uS) is illustrated in 
Fig.[TJa). Clearly, the equilibrium A^(to) reproduces cor- 
rectly the well known features of SIAM, such as the Hub- 
bard peaks at around uj = e and uj = e+U, and the Kondo 
peak centered at w = p cq = 0. 65 In Ref.|48|, the equilib- 
rium A^ui) of SIAM in the Kondo regime has been in- 
vestigated with the HEOM approach thoroughly and the 
existence of Kondo resonance is manifested by the correct 
universal scaling behavior there. In the nonequilibrium 
situation where an external voltage is applied antisym- 
metrically to the L and R reservoirs, the Hubbard peaks 
remain largely unchanged in both position and height. In 
contrast, the Kondo peak is split by the voltage into two, 
which appear at lo = A anc j w — _A anc [ correspond 
to the shifted reservoir chemical potentials /!l and /!r, 
respectively. Obviously, as the bias voltage A increases 
from to 0.2 meV, then to 0.7 meV, the progressive split- 
ting of Kondo peak is observed in Fig.QJa). Figure [ljb) 
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FIG. 1: HEOM calculated spectral function A(uj) = A r (u>) = 
A_i(lu) in unit of (-nT) -1 . The parameters adopted are the 
same as those specified in caption of Table |U (a) Illustration 
of the split process of equilibrium Kondo peak by the bias 
voltage A. (b) Temperature dependence of Kondo and Hub- 
bard peaks. The vanishment of peaks at the position of /il 
and /j,r manifests their Kondo nature. 



plots the calculated A(u>) of the same SIAM system at 
various temperatures. Apparently, as the temperature 
increases over an order of magnitude, the two Hubbard 
resonance peaks at uj = e and u> = e + U almost remain 
intact. In contrast, the split peaks at uj = fi^ and /i^ 
vanish quickly at the higher temperature. This clearly 
highlights the strong electron correlation features in the 
present noncquilibrium SIAM. 

The Hubbard peaks at uj = e and uj = e + U converge 
more rapidly than the resonance peaks at w = ±-y when 
truncation level L increases. This further confirms the 
Kondo nature of resonance peaks at uj = ±A6§. More- 
over, there is a long-time oscillatory tail in the real time 
evolution which is crucial for the appearance of Kondo 
peaks. We also find that the short-time dynamics of 
the retarded Green's function and high-frequency part 
of A(uj) converge more rapidly when the truncation level 
L increases. In other words, in the HEOM framework, 
one can extract the spectral function at high frequency 
range at a relatively lower truncation level and relatively 
shorter evolution time than those at resonance frequen- 
cies, maintaining quantitative accuracy. 

We then exemplify the numerical tractability of HEOM 
approach via evaluation of three types of response proper- 





1 1 1 1 




A=0 




A=0.7meV 







CO (meV) 

FIG. 2: HEOM calculated local charge fluctuation spectrum, 
Sq(uj) of Eq. 1)31 p. in unit of e 2 . The parameters adopted are 
the same as those specified in caption of Table U 

ties. These include the local charge fluctuation spectrum 
Sq(lu), local magnetic susceptibility xm{u), and differ- 
ential admittance spectrum G aa > (uj), defined respectively 
as follows, 

/oo 
dte iut ({AQ(t),AQ(0)}) st , (31) 
-OO 

poo 

Xm(u) = i / die*"* ( [M(t), M(0)] > st , (32) 
Jo 

POO 

G aa ,{uj) = -i dte iut ([i a (t),N a ,(0)]) (33) 
Jo 

In Eq. (EJ), AQ(t) = Q(t)-(Q) st , with Q = ^ n M being 
the total impurity occupation number operator. There- 
fore, AQ(t) describes the fluctuation of occupation num- 
ber around the averaged value. For xm(w) of Eq. (|3"2"|) . 
M = q^lbSz is the impurity magnetization operator, 
which originates from the electron spin polarization in- 
duced by external magnetic field. Here, g is the electron 
g-factor, [Ib is the Bohr magneton, and S z = {n^ — n±)/2 
is the impurity spin polarization operator. In Eq. (|33|) . 
G aa i (w) is just the half- Fourier transform of current- 
number response function of Eq. (f25|) or ((28)) . The time 
t = in the individual Eqs. (f3~T j) - (|3"3"|) represents the in- 
stant at which the external perturbation (magnetic field 
or bias voltage) is interrogated. 

It is worth pointing out that all the three types of 
response properties satisfy the following symmetry: the 
real (imaginary) part is an even (odd) function of oj. This 
is due to the time-reversal symmetry of steady-state cor- 
relation functions, i.e., Cab{^) = [Cba{— £)]*• In partic- 
ular, for Sq(w) of Eq. $5T$, A = B = AQ. Consequently, 
it can be shown that Sq(uj) is a real function. For clarity, 
Figs. [31 and [5] will only exhibit the uj > part of the 
dynamic response properties, while the u < part can 
be retrieved easily by applying the symmetry relation. 

The local charge fluctuation spectrum Sq(uj) has been 
studied in the context of shot noise of quantum dot 
systems ^2r— Figure [2] depicts the HEOM calculated 
Sq(ui) of the SIAM under our investigation. At equilib- 
rium, the spectrum exhibits a crossover behavior, where 
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FIG. 3: HEOM calculated local magnetic susceptibility, 
Xm(ui) of Eg. (|32|l . in unit of g 2 fi 2 B /kBT. The parameters 
adopted are the same as those specified in caption of Table [T] 



the two peaks centered at around uj = \e\ = 0.5 meV and 
u) = |e + U\ = 1 meV largely overlap and form a broad 
peak. The positions of these two peaks correspond to 
the excitation energies associated with change of impu- 
rity occupancy state. In nonequilibrium situation, the 
crossover peak is observed to move to a lower energy, 
since the chemical potential of reservoir R is drawn closer 
to the impurity state by the applied voltage. 

The local magnetic susceptibility \m{u) is a response 
property of fundamental significance, particularly for 
strongly correlated quantum impurity systems. It has 
been studied by various methods such as NRG.— Fig- 
ure |3] displays the HEOM calculated xa/M of the SIAM 
of our concern. In both equilibrium and nonequilibrium 
situations, the main characteristic features of xm(w) ap- 
pear at around zero energy. Moreover, the magnitude 
of XAf(w) is found to reduce significantly in presence of 
applied bias voltage, especially in the low energy range. 
This is consistent with the diminishing spectral density at 
around zero-frequency due to the voltage-induced split- 
ting of Kondo peak; see Fig.[TJ 

To verify the numerical accuracy of our calculated local 
magnetic susceptibility, we compare the HEOM approach 
with the latest high-level NRG method. The comparison 
is shown in Fig. 21 where the equilibrium static magnetic 
susceptibilities, xm{u = 0), of various symmetric SIAM 
systems are calculated to reproduce the Fig. 6 of Ref. [7(J 
Apparently, the HEOM and NRG results agree quantita- 
tively with each other. 

The dynamic admittance is one of the most exten- 
sively studied response properties of quantum dot sys- 
tems. The frequency-dependent admittance has been 
studied by scattering theory 7 -^— and nonequilibrium 
Green's function methodi 7 -^— The HEOM approach has 
also been used to calculate the dynamic admittance of 
noninteracting^i and interacting quantum dots . 52 ' 57 This 
was realized by calculating the transient current in re- 
sponse to a delta-pulse voltage^ In the following, we re- 
visit the evaluation for dynamic admittance G{u>) via an 
alternative route, i.e., by calculating the current-number 
response functions of Eq. (|25|) . 




T/T K 

FIG. 4: Static local magnetic susceptibility multiplied by tem- 
perature, Xm(<^> = 0) T (in unit of g 2 fj? B /kB), versus T/Tk 
for a series of equilibrium symmetric SIAM systems of dif- 
ferent U. Here, Tk is the Kondo temperature, and U, T, 
and W are in unit of F. The HEOM results (scattered sym- 
bols) are compared with the latest full density matrix NRG 
calculations (lines) of Ref. the Fig. 6 there. The NRG 
calculations are for very large reservoir bandwidth W, while 
the HEOM results are obtained with relatively smaller band- 
widths (W = 10 and W = 20) for saving computational cost. 
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FIG. 5: HEOM calculated dynamic admittance, G(uj) = 
j[Gll(lj) + Glr(u) -Gr-l(oj) -Gr-r.(oj)] with Gw(w) defined 
by Eq. (55) , in unit of e 2 /h. The parameters adopted are the 
same as those specified in caption of Table U 



Figure [5] depicts the HEOM calculated differential ad- 
mittance of the SIAM under study, G(w) = j[G\ J ^{ l -o) + 
Glr(w) - Gri» - GrrM] = \\G L (d) - G r (lj)}; cf. 
Eq. (JlOl), with G aa >(u) defined by Eq. (j33]). Here we 
have chosen antisymmetrically applied probe ac bias, 
<5Al(£) = —5An(t) = ^SA(t). As discussed extensively 
in Ref.|46|, the characteristic features of G(uj) appearing 
at around lu = |e| and uj = |e + XJ\ correspond to the 
transitions between Fock states of different occupancy, 
while the low-frequency features highlight the presence 
of dynamic Kondo transition. Apparently, the dynamic 
Kondo transition is suppressed by the applied voltage, 
which is analogous to the scenario of xm(w) as shown in 
Fig.El 



V. CONCLUDING REMARKS 

In this work, we develop a hierarchical dynamics ap- 
proach for evaluation of nonequilibrium dynamic re- 
sponse properties of quantum impurity systems. It is 
based on a hierarchical equations of motion formalism, 
in conjunction with a linear response theory established 
in the hierarchical Liouville space. It provides a unified 
approach for arbitrary response and correlation functions 
of local impurity systems, and transport current related 
response properties. 

The proposed hierarchical Liouville-space approach 
resolves nonperturbatively the combined effects of e- 
e interactions, impurity-reservoir dissipation, and non- 
Mar kovian memory of reservoirs. It provides a uni- 
fied formalism for equilibrium and nonequilibrium dy- 
namic response properties of quantum impurity systems 
and can be applied to more complex quantum impu- 
rity systems without extra derivation efforts. Moreover, 
the HEOM results converge rapidly and uniformly with 
higher-tier ADOs included explicitly and one can often 
obtain convergent and quantitatively accurate results at 
a relative low truncation level L. For equilibrium proper- 
ties, the same level of accuracy as the latest state-of-the- 
art NRG method can be achieved via HEOM approach. 

In conclusion, the developed hierarchical Liouville- 
space approach provides an accurate and universal tool 
for investigation of general dynamic response properties 
of quantum impurity systems. In particular, it addresses 
the nonequilibrium situations and resolves the full fre- 
quency dependence details at quantitative accuracy. It 
is thus potentially useful for exploration of quantum im- 
purity systems and strongly correlated lattice systems 
(combined with dynamical mean field theory), which are 
of experimental significance for the advancement of na- 
noelectronics, spintronics, and quantum information and 
computation. 
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Appendix A: Supplement Materials 

In this Supplemental Material, we will discuss some 
details of the numerical implementation of the proposed 
HEOM formalism to the nonequilibrium dynamical prop- 
erties of quantum impurity systems. 

(i) In the main text, the spectral function of a- 
reservoir assumes a diagonal and Lorentzian form, i.e., 
r 

Ja^M = Spv ( ^_ M ° )2 ° , with r Q and W a being 
the linewidth and bandwidth parameters, respectively. 
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L, M 




P-n 


I (nA) 


1, 9 


0.500 {0.500) 


0.001 (0.000) 


0.003 


2, 9 


0.441 (0.462) 


0.025 (0.027) 


4.654 


3, 9 


0.439 ( 0.454) 


0.024 (0.025) 


4.920 


4, 9 


0.440 (0.457) 


0.024 ( 0.024) 


4.799 


4, 10 


0.440 (0.457) 


0.0237 (0.024) 


4.805 


4, 11 


0.440 (0.457) 


0.0237 (0.0238) 


4.806 


5, 9 


0.440 (0.457) 


0.024 ( 0.024) 


4.799 


5, 10 


0.440 (0.457) 


0.0237 (0.024) 


4.805 


5, 11 


0.440 (0.457) 


0.0237 (0.0238) 


4.806 



TABLE II: Spin-/i single- and double-occupation probabili- 
ties (Pfj, and P-C4.), and steady-state current of an SI AM with 
two electrons reservoirs under an antisymmetrically applied 
bias voltage of A = —V = 0.7 meV. Calculations are done 
by solving the HEOM truncated at different level L with dif- 
ferent number of exponential terms M. The parameters are 
the same as those specified in caption of Table I of the main 
text. For comparison, the numbers of equilibrium situation 
of V — are shown in the parentheses. 
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CO (meV) 

FIG. 6: The convergence behavior of HEOM calculated spec- 
tral function A(uj) of SIAM (in unit of I/71T) with respect to 
the number of exponential terms M in spectrum decomposi- 
tion. Here the truncation level is fixed at L = 4. The param- 
eters are the same as those specified in caption of Table I of 
the main text. 



Then in the contour integration evaluation of it) = 
r fa e Tl J hi'} U uT\ a total number of M = N + 1 

j — 00 l-\-e p k^ — f-a ) 1 

poles will be involved, including one Drude pole of the 
reservoir spectral density and N poles of the Fermi dis- 
tribution function. 

In general, the number K of distinct ADO indexes in 
the HEOM proposed in the main text amounts to K — 
2N a N fJi M, with being the number spin-orbitals of 
system in direct contact to leads. The factor 2 accounts 
for the two choices of the sign variable tr, while N a = 2 
for the distinct a = L and R leads. 

We have known the fact that overall computational 
cost increases dramatically with both L and K. To mini- 
mize the computational cost while maintaining the quan- 
titative accuracy, a highly efficient reservoir spectrum 
decomposition scheme is needed to optimize the size of 
K . Various decomposition schemes have been developed, 
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including Matsubara spectrum decomposition^ partial 
fractional decomposition^ Pade spectrum decomposi- 
tion (PSD),^& and hybrid schemed To our knowledge, 
PSD scheme has the best performance and acquires the 
same accuracy with the minimal size of K until now. In 
the present work, the [N — l/N] Pade spectrum decom- 
position schem e 55 ! 56 is used for the efficient construction 
of the hierarchical Liouville propagator. 

There is an overlap between Fermi distribution func- 
tion fa — Ma) an d bath spectral density J^av ( w — Ma ) 
in the definition of bath correlation function C£jf£(t) = 

f^due^f^iu - (iaiJZ^iu - Ma). Due to the fi- 
nite bandwidth of reservoirs in realistic quantum impu- 
rity system, the residual between f!^\u> — Ma) an d its 
Pade series expansion (with large enough N) at u » W a 
will not introduce any approximation. Therefore, the ex- 
ponential series expansion of bath correlation function 
with large enough M can efficiently resolves the non- 
Markovian memory of reservoirs in both equilibrium and 
nonequilibrium situations. 

To obtain quantitatively accurate numerical results, we 
increase the number of exponential terms M in spectrum 
decomposition continually until convergent results are ar- 
rived in practiced 

Table |TT] lists the probabilities that the impurity is 
singly occupied by spin-// electron (P M = (p\p(T,V)\p) 
with fx =f or I); or doubly occupied (P^ = 
{U\p(T,V)\U))- Here, p(T,V) is the nonequilibrium 
steady-state reduced density matrix under temperature 
T and antisymmetric applied voltage V. Calculations are 
done at different truncation level L and different number 
of exponential terms M (up to L = 5 and M = 11). Ap- 
parently, the HEOM results converge rapidly with the 
increasing L, i.e., with higher-tier ADOs included ex- 
plicitly in HEOM. In particular, the remaining relative 
deviations between the results of L = 4 and L = 5 are less 
than 0.1%, indicating that the L = 4 level of truncation is 
sufficient for the present set of parameters. At the same 
time, one can observe that M = 9 is sufficient to yield 
convergent results here. These are further affirmed by 
the calculated steady-state current /(V) across the im- 
purity, which also converges quantitatively with rather 
minor residual uncertainty at L = 4 and M = 9. 

Fig-El further illustrates the convergence behavior of 
A(us) of SI AM with respect to the number of exponential 
terms M at the fixed truncation level L = 4. Obviously 
71/ = 9 is large enough for the adopted parameters of the 
numerical examples in the main text. 

(ii) The terminal truncation level L max is often too 
high to reach in practical implementation, and a trun- 
cation at a relatively low level L is inevitable. Here 
we adopt a straightforward truncation scheme, i.e., set 
all the higher-tier ADOs (n > L) zero directly. And 
we increase the truncation level L continually until the 
convergent numerical results are obtained. For quantum 
impurity systems with nonzero e-e interactions, calcula- 
tions often converge rapidly and uniformly with the in- 
creasing truncation level L. For weak to medium system- 
reservoir coupling strength, one can often obtain quan- 
titatively accurate results at a relatively low L. HEOM 
calculated spectral function A(u) of SIAM at different 




(£> (meV) 

FIG. 7: The convergence behavior of HEOM calculated spec- 
tral function A(us) of SIAM (in unit of I/71T) with respect to 
the truncation level L = 2,3,4 with the fixed number of expo- 
nential terms M = 9. The parameters are the same as those 
specified in caption of Table I of the main text. 




t(ps) 

FIG. 8: The convergence behavior of HEOM calculated 
G£*(i > 0) of SIAM with respect to the truncation level 
L = 2,3, 4, with fixed M = 9. The parameters are the same 
as those specified in caption of Table I of the main text. The 
inset magnifies the long time oscillatory tail in real time evo- 
lution. 



truncation level L = 2, 3, 4 with the fixed number of ex- 
ponential terms in spectrum decomposition (M = 9) are 
plotted together in Fig. [3 to illustrate the convergence 
behavior and gain the feeling on the numerical aspects 
of the proposed formalism. Obviously, the numerical re- 
sults converge rapidly and uniformly with the increasing 
truncation level L. Fig. [7] further confirms the fact that 
L = 4 level of truncation is sufficient for the set of pa- 
rameters adopted in the main text. One can also observe 
that the resonance Hubbard peaks at around cj = e and 
u = e + U converge more easily than the split Kondo 
peaks at oj = ±y with increasing truncation level L. 
This further confirms the strong correlation nature of 
resonance at us = ±-y- Moreover, the numerical results 
at high and far-from-resonance frequency range converge 
more easily than those at Kondo resonance frequencies. 

(iii) The hierarchical Liouvillc-space linear response 
theory proposed in the main text can be implemented 
in either time-domain or frequency domain. In the time- 
domain scheme, one first carry out the real-time evo- 
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lution under the unperturbed HEOM propagator Q s (t), 
starting from the initial condition /5(0), i.e., 

p(t) = g s (t)p(o)- (Ai) 

Then the correlation/response functions can be extracted 
from certain ADOs components {pj^ j (t)} in p(i), as 
in Eq. (20) of the main text. At last the frequency- 
dependent dynamical properties are obtained straight- 
forwardly by a half Fourier transform. 

Fig. [8] depicts the HEOM calculated real part of 
GW(t > 0) = -i({a^(t),a\(0)}} st of SIAM system at 
different truncation level L = 2, 3, 4. The number of 
exponential terms in spectrum decomposition is fixed at 
M = 9. The imaginary part of G^(t > 0) has the same 
convergence behavior and is ignored here. 

Obviously the short-time dynamics converges more 
easily than long-time dynamics when the truncation level 
L increases. This is consistent with the frequency-domain 
convergence features with respect to truncation level in- 
dicated in Fig. [7] This is because that dynamical proper- 
ties at high and far-from-resonance frequency range are 
dominated by the short-time evolution. 

As illustrated in the inset of Fig. HI the real time evolu- 
tion has a long-time oscillatory tail. To obtain quantita- 
tive accurate spectral function at high and off-resonant 
frequencies, one only need evolve relative short time t. 
In contrast, the accurate Kondo peaks can be achieved 
only when long enough oscillatory evolution is included. 
Therefore the long-time oscillatory tail is crucial for the 
Kondo resonance peaks. It can be viewed as another 
evidence of the existence of Kondo correlation. 

Therefore for strongly correlated quantum impurity 
systems at low temperatures, where Kondo signatures ap- 
pear at the position of chemical potentials p, a , one must 
propagate G^(t > 0) to a long enough time to include 
the significant long-time memory of reservoir. Conse- 
quently, it is very time-consuming to obtain quantita- 
tive accuracy at the frequencies of Kondo peaks. For the 
present parameters in the main text, when A = 0.7meV, 
L = 4 and AI = 9, one need to evolve about 190ps (con- 
suming about 10000 minutes CPU time) to obtain the 
accurate Kondo peaks. 

An alternative numerical scheme is carried out in fre- 
quency domain. In frequency-domain scheme one solves 
the following sparse linear algebra problem at a fixed fre- 
quency u>, 

n 

- ,„ H - *. (0) = - [*C + E iir\ H 

n 

j r—1 

(A2) 

where Pj^:..j (w) is the half- Fourier transformation of 



p\\...j n {t)- And p\\...j n (0) is the initial condition of the 
corresponding correlation/response function. To obtain 
the whole spectrum, a full scan for the whole frequency 
domain is needed. 

The equivalence between the time- and frequency- 
domain schemes is exemplified with Fig.[§l where the 
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FIG. 9: HEOM calculated spectral function A(u) of SIAM 
obtained in both time domain with t — 190ps and frequency 
domain, in unit of (71T) -1 . The parameters are the same as 
those specified in Fig.[8] The truncation level is fixed at L — 4 
and the number of exponential terms in spectrum decompo- 
sition is fixed at M = 9. 



spectral function A(u) of a SIAM is calculated with both 
time- and frequency-domain schemes. Apparently, the 
results of both schemes agree perfectly with each other. 

In practice, one can employ a hybrid time- and 
frequency-domain scheme. The overall profile can be 
plotted with the time-domain scheme, while the detailed 
structures of the resonance and Kondo peaks can be ex- 
ploited with the frequency-domain scheme. This hybrid 
scheme provides an accurate and efficient method to cal- 
culate the full spectrum of a strongly correlated quantum 
impurity system in the Kondo regime. 

In summary, the truncation level and the number of 
exponential terms in spectrum decomposition are contin- 
uously increased until convergent numerical results are 
obtained. Therefore the proposed HEOM nonequilib- 
rium linear response theory provides an quantitatively 
accurate numerical tool for arbitrary response and corre- 
lation functions of local impurity systems, and transport 
related response properties. In the frame of our HEOM 
formalism, one can extract dynamical properties at high 
and off-resonant frequency range at relative lower trunca- 
tion level and relative shorter evolution time than those 
at resonance frequencies, maintaining quantitative accu- 
racy. The hybrid time- and frequency-domain scheme 
provides an accurate and efficient method to obtain the 
full spectrum of a quantum impurity system in the Kondo 
regime. 
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